function pso_optimization()
    % PSO参数设置
    n_particles = 100;     % 粒子数量
    max_iter = 30;         % 最大迭代次数
    w = 0.7;               % 惯性权重
    c1 = 1.5;              % 个体学习因子
    c2 = 1.5;              % 群体学习因子
    
    % 变量边界(与GA相同)
    lb = [0.10; 0.05; 0.01; 0.001; 0.05];
    ub = [0.11; 0.3; 0.1; 0.05; 0.2];
    dim = length(lb);      % 变量维度
    
    % 初始化粒子群
    particles = zeros(n_particles, dim);
    velocities = zeros(n_particles, dim);
    pbest_pos = zeros(n_particles, dim);
    pbest_val = inf(n_particles, 1);
    gbest_pos = zeros(1, dim);
    gbest_val = inf;
    
    % 记录收敛过程
    convergence = zeros(max_iter, 1);
    
    % 初始化粒子位置和速度
    for i = 1:n_particles
        particles(i,:) = lb' + (ub' - lb') .* rand(1, dim);
        velocities(i,:) = 0.1 * (ub' - lb') .* randn(1, dim);
        pbest_pos(i,:) = particles(i,:);
        
        % 计算初始适应度
        current_val = hmm(particles(i,:));
        pbest_val(i) = current_val;
        
        % 更新全局最优
        if current_val < gbest_val
            gbest_val = current_val;
            gbest_pos = particles(i,:);
        end
    end
    
    % PSO主循环
    for iter = 1:max_iter
        for i = 1:n_particles
            % 更新速度
            r1 = rand(1, dim);
            r2 = rand(1, dim);
            velocities(i,:) = w * velocities(i,:) + ...
                             c1 * r1 .* (pbest_pos(i,:) - particles(i,:)) + ...
                             c2 * r2 .* (gbest_pos - particles(i,:));
            
            % 限制速度范围
            max_velocity = 0.2 * (ub' - lb');  % 确保是行向量
            % 逐个元素限制速度
            for d = 1:dim
                if velocities(i,d) > max_velocity(d)
                    velocities(i,d) = max_velocity(d);
                elseif velocities(i,d) < -max_velocity(d)
                    velocities(i,d) = -max_velocity(d);
                end
            end
            
            % 更新位置
            particles(i,:) = particles(i,:) + velocities(i,:);
            
            % 确保位置在边界内
            particles(i,:) = min(max(particles(i,:), lb'), ub');
            
            % 计算适应度
            current_val = hmm(particles(i,:));
            
            % 更新个体最优
            if current_val < pbest_val(i)
                pbest_val(i) = current_val;
                pbest_pos(i,:) = particles(i,:);
                
                % 更新全局最优
                if current_val < gbest_val
                    gbest_val = current_val;
                    gbest_pos = particles(i,:);
                end
            end
        end
        
        % 记录当前最优值
        convergence(iter) = gbest_val;
        
        % 显示迭代信息
        fprintf('Iteration %d: Best Value = %.4f\n', iter, gbest_val);
        
        % 绘制收敛曲线
        figure(1);
        plot(1:iter, -convergence(1:iter), 'b-', 'LineWidth', 2);
        xlabel('Iteration');
        ylabel('Best Absorption');
        title('PSO Convergence');
        grid on;
        drawnow;
    end
    
    % 显示最终结果
    fprintf('\nOptimization completed!\n');
    fprintf('Best thickness found:\n');
    fprintf('t1 = %.4f um\n', gbest_pos(1));
    fprintf('t2 = %.4f um\n', gbest_pos(2));
    fprintf('t3 = %.4f um\n', gbest_pos(3));
    fprintf('t4 = %.4f um\n', gbest_pos(4));
    fprintf('t5 = %.4f um\n', gbest_pos(5));
    fprintf('Best absorption = %.2f%%\n', -gbest_val*100);
    
    % 保存结果
    save('pso_results.mat', 'gbest_pos', 'gbest_val', 'convergence');
end